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Introduction 


The  analysis  of  large  structures,  especially  in  three  dimensions,  can  result  in 
matrices  that  demand  an  exceptionally  large  amount  of  computer  storage.  The  storage  needs  of 
these  matrices  depend  to  a  large  extent  on  their  sparsity  and  the  data  structure  that  is  used  u. 
store  them.  The  choice  of  the  data  structure  in  turn  depends  on  the  method  that  is  used  ;  j  so. 
the  associated  system  of  equations.  Presently,  most  solution  schemes  used  in  finite  element  com¬ 
puter  programs  are  based  on  direct  methods,  i.e.  triangular  factorisation  of  the  stiffness  matrix, 
K.  Starting  from  a  given  mesh  description,  a  finite  element  program  performs  the  following  steps: 

1.  determine  the  sparsity  structure  of  K, 

2.  renumber  the  equations  to  reduce  the  storage  demands  of  K, 

3.  reserve  the  required  storage  for  K, 

4.  compute  the  element  matrices, 

5.  assemble  the  element  matrices  into  K, 

6.  compute  the  triangular  factorization  of  K, 

7.  solve  the  associate  system  of  equations. 

In  many  applications  the  available  primary  memory  is  not  sufficient  to  store  the  assembled 
matrix,  and  therefore  secondary  storage  is  used.  In  this  circumstance,  steps  4  through  7  involve 
data  transfers  between  primary  and  secondary  store,  often  referred  to  as  I/O.  In  this  case,  K  is 
partitioned  into  blocks  and  each  block  is  assembled  and  stored  on  secondary  store.  The  blocks 
are  then  brought  back  into  the  main  memory  to  form  the  factors  of  K.  For  large  enough  prob¬ 
lems  the  I/O  costs  can  dominate  the  computation  costs. 

A  great  deal  of  effort  has  been  expended  to  develop  new  procedures  for  reordering  the  equili¬ 
brium  equations,  thus  reducing  the  overall  storage  requirements  in  the  solution  steps  5,  6  and  7 
[1,2,3).  This  is  motivated  by  the  fact  that  reduction  in  storage  translates  directly  into  a  savings 
in  the  I/O  costs.  Among  the  many  solution  schemes  used,  the  frontal  method  |4j  and  the  profile 
or  skyline  method  (1|  are  probably  the  most  popular.  In  [8|  it  is  shown  that  when  the  same  nodal 


elimination  ordering  is  used  the  profile  method  performs  the  same  number  of  operations  a*  tor 
frontal  method;  in  [5]  an  algorithm  is  described  that  delivers  a  good  frontal  node  ordering  for  the 
profile  method.  The  significant  difference  between  these  methods  is  that  the  frontal  method  olteii 
combines  steps  A,  5  and  6.  In  this  way  the  I/O  during  the  assembly  step  is  overlapped  with  the 
I/O  in  the  factorization  step;  thus,  the  frontal  method  results  in  a  saving  that  is  equal  to  the  co» 
of  I/O  in  step  4.  Alternatively,  the  assembly  of  a  profile  stored  matrix  which  is  partitioned  into 
blocks  requires  a  multiple  pass  through  the  elements  to  perform  the  assembly.  The  principle 
differences  in  the  I/O  costs  of  the  two  methods  during  factorization  may  be  traced  to  the  above 
differences  in  the  partitioning  of  the  matrix  into  a  frontal  or  a  profile  form.  The  principle  disad 
vantage  of  the  frontal  method  is  the  added  overhead  to  retain  a  small  front  width  during  the  tri¬ 
angular  factorization  in  step  6  and  subsequent  resolutions  in  step  7.  For  example,  in  resolution  of 
equations  this  added  overhead  may  lead  to  CPU  C06ts  which  are  several  times  those  of  a  profile 
stored  resolution. 

In  this  paper  we  take  a  different  approach.  We  use  the  following  simple  observation: 

The  assembly  process  la  Independent  of  the  solution  procedure. 

In  other  words,  one  should  use  the  most  efficient  data  structure  for  the  assembly  process,  step  5, 
and  then  restructure  the  data  for  ones  favorite  solution  scheme,  i.e.,  either  frontal  or  profile.  In 
this  way  one  can  achieve  the  same  reductions  in  I/O  as  the  frontal  method  and  at  the  same  time 
maintain  high  modularity  of  the  program.  Here  we  develop  a  data  structure  that  stores  only  the 
nonzero  terms  in  the  stiffness  matrix  in  a  compacted  form,  and  present  an  algorithm  for  the 
assembly  of  K  for  this  storage  method.  This  approach  results  in  considerable  reduction  in  the 
storage  needs  during  the  assembly  process.  Therefore  large  matrices  often  may  be  fully  assembled 
in-core  resulting  in  a  considerable  reduction  in  I/O. 

This  approach  has  the  added  advantage  that  the  program  is  not  built  around  a  single  equa¬ 
tion  solver.  One  can  have  many  solution  procedures  by  simply  expanding  the  compacted  struc¬ 
ture  of  K  into  a  form  appropriate  for  each  particular  solution  method.  Furthermore,  the  com¬ 
pacted  structure  can  be  used  directly  for  iterative  solution  techniques  such  as  the  conjugate 


gradient  type  methods  |6,7]. 


Storage  Scheme 

We  now  describe  the  compacted  structure  that  is  used  to  store  K.  We  only  consider  sym¬ 
metric  matrices,  although  the  extension  to  the  nonsymmetrk  case  is  trivial,  and  store  only  u.- 
upper  triangular  part  of  K  is  stored.  The  diagonal  terms  of  K  will  be  stored  separately  in  a  *ia 
gle  array  of  length  n,  where  n  is  the  total  number  of  equations.  The  remaining  off-diagonals  will 
be  placed  in  a  second  array  of  length  r,  where  r  is  the  total  number  of  nonzero  off-diagonal  terms 
in  the  upper  triangular  part  of  K.  All  the  elements  in  the  same  column  will  be  placed  consecu¬ 
tively  in  this  array,  starting  from  the  top  of  the  column  down  to  the  diagonal  (excluding  the  diag¬ 
onal  term).  The  columns  are  stored  consecutively  from  the  second  to  the  last.  For  each  entry  in 
this  array  we  store  its  row  number  in  a  corresponding  integer  array  of  length  r.  The  example 
below  demonstrates  the  final  storage  scheme. 

Example: 

Consider  the  matrix 


1  9 

2  10  11 
3 

4 


12 

15 

13 

14 

16 

17 

19 

21 

5 

18 

20 

6 

22 

7 

23 

8 


The  array  Hag  contains  the  diagonal  terms  of  K  as  shown  below. 


Hag  -  11,2,3,4,5,6,7,8)  . 

A  real  array  then  stores  the  off-diagonal  terms  of  K  and  a  corresponding  integer  array  denotes  the 
row  number  of  each  off-diagonal  term  as  show  below. 

off-iiat  -  (9,10, 11, 12, 16, 13, 17, 18, 14, 19, 20, 15,21,22,23| 
trow— (1,  2,  2,  1,  3,  2,  3,  5,  2,  4,  5,  1,  4,  6,  7) 

la  addition,  a  single  integer  array  of  length  n  is  also  required  to  point  to  the  end  of  entries  from  a  j 

given  column.  For  the  above  example  this  army  is  ] 
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jcol  «  (0,1,2,3,5,8,11,15) 

The  total  storage  requirement  is  r  +  n  real  words  and  r  +  n  integer  words.  Using  a  16 

bit  integer  word,  ^-(n  +  r)  real  words  (64  bit)  will  be  sufficient.  Then  the  largest  number  of 
4 

equations  that  can  be  solved  this  way  is  216  -  1  ~  32000.  With  a  32  bit  integer  word,  the 
3 

storage  required  will  be  —(n  +  r)  real  words. 

Derivation  of  the  Assembly  Process 

In  this  section  we  give  a  step  by  step  derivation  of  the  assembly  algorithm.  Each  step  is 
demonstrated  with  the  aid  of  the  mesh  example  in  Pig.  1.  First  we  introduce  some  notation.  A 
finite  element  mesh  is  denoted  by  M  =  [E  ,  N)  where  E  and  N  represent  the  collection  of  ele¬ 
ments  and  nodes  in  the  mesh.  A  part  of  the  input  information  provided  to  a  finite  element  pro¬ 
gram  is  the  set  of  nodes  belonging  to  an  element.  This  we  denote  as  N,QN.  For  example  ele¬ 
ment  4  in  Fig.  1  has  the  connectivity  set  N4  =  {5, 8, 0,6 ).  In  Table  1  we  give  the  complete  list 
of  the  connectivity  sets  Nt  for  each  element  in  example  1.  These  data  are  usually  assembled  in  a 
single  array  known  as  the  connectivity  array.  The  complete  set  {  Ne  V  e€E }  is  sufficient  to 
describe  the  connectivity  of  a  given  mesh.  Another  part  of  the  input  data  is  the  boundary  condi¬ 
tions  that  determine  the  set  of  the  indices  of  all  the  active  degrees  of  freedom  at  node  p.  We 
denote  this  set  as  Ur.  In  Table  1  we  give  the  set  (  Uf  Vp€Ar)  for  the  example  in  Fig.  1.  In  this 
example,  we  assume  that  there  are  two  degree  of  freedom  at  each  node.  The  collection  of  Nt, 
column  2,  and  Vr,  column  4,  given  in  Table  1  is  sufficient  to  determine  the  sparsity  structure  of 
the  stiffness  matrix  associated  with  a  given  mesh. 
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Element  Connectivity 

Boundary  Conditions 

Element 

Set  of  Nodes  for 

Node 

Set  of  Unknowns  for 

No.  e 

each  Element,  N, 

No.  p 

the  Active  nodes,  Vt 

1 

{  1.  4,  5,  2  } 

2 

{1} 

{  2,  5,  6,  3  } 

4 

{2} 

(  4,  7,  8,  5  } 

5 

{3,4  ) 

{  5,  8,  9,  6  } 

6 

m 

5 

{  7,  II,  8  } 

7 

{0} 

6 

{  10,  II,  7  } 

8 

{7,8} 

7 

{  8,  11,  9  } 

9 

{9} 

8 

{  11.  12,  9  } 

11 

{ io) 

Table  1.  Connectivity  sets  and  active  degrees  of  freedom  for  the  mesh  in  Fig.  1. 

Our  objective  here  is  to  find  the  set  of  indices  of  the  unknowns  that  are  coupled  with  a 
given  degree  of  freedom.  This  is  precisely  the  row  number  of  each  nonzero  term  in  a  given 
column  of  the  stiffness  matrix,  trow,  that  is  required  for  the  storage  scheme  described  in  the  pre¬ 
vious  section. 

First,  we  must  establish  the  set  of  elements  that  are  connected  to  each  node.  This  can  be 
done  by  inspecting  the  element  connectivity  sets.  Looking  at  the  second  column  of  Table  1,  for 
example  node  4  appears  twice,  in  rows  1  and  3.  We  then  conclude  that  node  4  is  connected  to 
elements  1  and  3.  This  process  must  be  repeated  for  each  node.  The  difficulty  here  is  that  we  do 
not  know  apriori  the  number  of  storage  locations  needed  to  identify  the  set  of  elements  for  each 
node.  For  this  reason  the  above  process  is  carried  out  in  two  steps.  The  number  of  elements  con¬ 
nected  to  each  node  is  determined  and  stored  first.  We  refer  to  this  as  the  E-degree  (element 
degree)  of  each  node.  In  the  example  the  E-degree  of  node  4  is  2.  The  E-degree  also  determines 
the  length  of  the  array  that  is  required  to  keep  the  set  of  elements  connected  to  each  node.  For 
each  node  p  we  denote  this  set  by  Ef  C  E.  Then  the  above  process  is  simply  to  evaluate 


Et  =  {t\pZN,)  (1) 

for  each  node  p.  This  equation  may  be  thought  of  as  finding  the  pseudo-inverse  of  the  connec¬ 
tivity  array-  The  E-degree  of  node  p  is  the  number  of  terms  in  Ef  (the  cardinality  of  Ef ),  See 


columns  two  and  three  of  Table  2  for  the  E-degree  and  the  complete  set  of  Ep  for  the  nodes  in 
example  1. 


Next,  we  find  for  the  set  of  nodes  that  are  adjacent  to  each  node  p .  We  denote  this  tj 


A, CN.  This  is  the  set  of  all  nodes  that  belong  to  an  element  with  p  as  one  of  its  nodes  H  <<nr  r. 


established  the  set  of  elements  connected  to  node  p,  the  adjacent  nodes  are  all  the  ot!  . . 


belonging  to  these  elements.  In  example  1  node  4  is  connected  to  elements  1  and  3  (see  column  2 


of  Table  2).  The  set  of  nodes  belonging  to  elements  1  and  3  are  obtained  by  inspecting  column  2 


of  Table  1;  these  are  {  1,  4,  S,  2  }  and  {  4,  7,  8,  5  }  respectively.  Then  the  set  of  nodes  adjacent 
to  4  is  —  {  1,  5,  2,  7,  8  }.  The  N-degree  (nodal  degree)  of  a  node  is  the  number  of  nodes  adja¬ 


cent  to  it  and  is  the  cardinality  of  Ar.  In  columns  4  and  5  of  Table  2  we  give  the  N-degree  and 
the  adjacency  set  of  each  node  in  the  nodal  graph  of  example  1  (Fig.  2).  This  step  is  simply  to 


evaluate  the  equation 


X,=  (J  Nt  -  n  (2) 

Note  that  both  Ap  and  N-degree  can  be  obtained  in  the  same  loop. 


Node 

P 

E 

degree 

E, 

N 

degree 

S 

degree 

1 

1 

{1} 

3 

{4,5,2} 

<t> 

{1.2} 

5 

{  1,  4,  5,  6,  3  } 

(1) 

{2} 

3 

{  2,  5,  6  } 

{2} 

{  1.  3  } 

5 

{  1,  2,  5,  7,  8  } 

{1.2} 

5 

{  1.  3,  4,  2  } 

8 

{  1,  4,  7,  8,  9,  6,  3,  2  } 

{  1,  4,  3,  2  } 

6 

(2,4) 

5 

{  2,  5,  8,  9,  3  } 

{  2,  5,  3  } 

7 

{  3,  5,  0  } 

5 

{  10,  11,  8,  5,  4  } 

{4,5} 

8 

{  3,  5,  7,  4  } 

6 

{  11,  9,  6,  5,4,7  } 

{  6,  5,  4,  7  } 

9 

3 

{4,7,8} 

5 

{  11,  12,  6,5,8  } 

{  6,  5,  8  } 

10 

1 

{«> 

2 

{  7,  11  } 

1 

{7} 

11 

4 

{  6,  5,  7,  8  } 

5 

{  10,  12,  9,  8,  7  } 

4 

{  10.  9,  8,  7  } 

12 

1 

2 

{».»} 

2 

{  11.9) 

Table  2. 


The  result  of  algorithm  for  establishing  the  row  index  of  the  nontero  terms  in  K 
for  example  1. 
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Since  we  want  to  store  only  the  upper  triangular  part  of  K  we  need  to  store  only  a  subset  of 
Af.  This  will  be  the  set  of  nodes  in  A,  with  an  index  less  than  p ;  that  is 
S,  =  {i  |  »€A,  and »  <  p  }.  We  refer  to  the  number  of  terms  in  St  as  the  S- degree  (  .  rr.'- 
degree)  of  a  node.  The  set  Sf  is  only  useful  when  the  numbering  of  the  unknowns  are  surh  thru 
when  i  <  j  all  the  unknowns  at  node  »  have  a  smaller  index  than  the  unkuowns  at  nMc 
Whenever  this  is  not  true  it  is  necessary  to  use  the  complete  set  of  adjacent  nodes  A,  together 
with  the  numbers  of  the  unknowns  for  each  node  (e.g.,  see  listing  in  Appendix  A.). 

Finally,  for  a  given  unknown  at  node  p  with  index  j£Uf  we  find  the  set  of  the  indices  of  all 
other  unknowns  that  are  coupled  with  u}.  This  will  be  the  set  of  row  indices  Rt  for  nonreros  in 
jT-th  column  of  K.  Then 

-  u  u>  (3) 

{•€>«,) 

For  example  1  E,  is  the  adjacency  set  of  j  in  the  graph  of  the  unknowns  in  Fig.  3.  Since  we  only 
store  the  upper  triangular  part  of  K  we  scan  through  S,  and  use  the  subset  defined  by: 

—  »€#,  and  »'<;'}  (4) 

R,  is  the  row  index  of  all  the  nonzero  terms  in  the  j-th  column  of  the  upper  triangular  part  of  K. 
The  complete  set  of  R ,  for  the  example  problem  is  presented  in  Table  3. 


D.O.F. 

; 

Node 

P 

l 

2 

{«} 

2 

4 

{U 

3 

5 

{1.2} 

4 

5 

{  3,  2,  1  } 

5 

6 

{  1  3,  1  } 

6 

7 

{  3,  4,  2  } 

7 

8 

{  6,  2,  4,  5,  3  } 

8 

8 

{  2,  3,  4,  5,  8,  7  } 

9 

9 

{  7,  4,  5,  3,  8  } 

10 

11 

{  9,  6,  7,  8  } 

Table  3.  The  row  indices  of  the  noniero  terms  in  the  upper  triangular  part  of  K. 


-  8- 


The  listing  of  n  FORTRAN  program  that  performs  all  tbe  steps  that  is  described  m  ibis 
section  is  given  in  Appendix  A.  In  this  Appendix  we  also  provide  tbe  subroutine  that  use-  the 
row  indices  to  perform  the  assembly  of  K. 

Numerical  Result 

We  use  tbe  algorithm  described  in  tbe  previous  section  to  assemble  the  stiffness  matrix  s  . 
tbe  problems  described  in  Table  6.  The  total  storage  required  during  the  assembly  step  is 
evaluated.  We  compare  these  results  to  similar  results  obtained  when  the  assembly  is  performed 
directly  into  a  profile  data  structure.  The  storage  requirement  of  the  compacted  assembly  is  not 
effected  by  the  node  ordering.  For  the  assembly  into  a  profile  form,  we  numbered  the  nodes 
across  the  width  of  the  mesh  to  reduce  the  bandwidth  of  the  stiffness  matrix.  Although,  the 
bandwidth  could  have  been  reduced  further  using  a  renumbering  scheme  such  as  [2.3],  we  omitted 
this  step  for  simplicity.  The  results  for  examples  2  to  6,  given  in  Table  6,  are  presented  in  Table 
4  below. 


Description 
of  Problem 

No.  of 
nodes 

No.  of 
elements 

No.  of 
equations 

Stores  for 
profile 

Stores  for 
Compacted  K 

Cantilever  type  structure 

225 

184 

428 

10204 

5340 

Small  Cylinder  structure 

231 

200 

440 

10492 

5679 

Large  Cylinder  structure 

496 

450 

960 

32542 

12715 

4X4X4  solid  structure 

125 

64 

I  '  ' 

21945 

11634 

8X8X8  solid  structure 

729 

512 

1944 

470043 

94272 

Table  4.  Comparison  of  the  storage  demands  of  profile  and  compacted  assembly  for  K  in 
examples  2  to  0. 

The  results  in  Table  4  is  obtained  based  on  the  assumption  that  a  real  word  is  twice  as  long  as  an 
integer  word.  We  observe  a  reduction  from  40%  for  two  dimensional  (2-D)  problems  to  80%  for 
3>D  problems  for  these  examples.  The  reductions  will  be  more  if  short  integer  words  are  used.  It 
is  interesting  to  note  that  the  required  storage  for  compacted  structure  varies  linearly  with  the 
number  of  equations.  Therefore,  the  saving  will  be  more  for  larger  problems.  In  Table  5  we  give 


-li¬ 


the  storage  counts  for  the  two  methods  considered  here  on  square  mesh  in  2-1)  and  cnl  cs  :u 
a  function  of  the  number  of  equations.  To  obtain  these  estimates  we  assumed  that  there  is  >j»fy 
one  degree  of  freedom  per  node. 


Dimension 
of  Problem 

Half 

bandwidth 

Profile 

storage 

Compacted 

storage 

1 

2 

5/2  n 

3n 

2 

„«/* 

ns* 

15/2  n 

3 

21  n 

Table  5.  Estimated  storage  needs  for  each  scheme  on  regular  mesh,  n  is  the  number  of  equations. 
Conclusion 

The  essential  steps  in  a  finite  element  program  can  be  modified  to  make  use  of  the  com* 
p acted  assembly  described  here.  Accordingly,  we  perform  the  following  steps: 

1.  Obtain  the  row  indices  of  the  nontero  terms  in  K. 

2.  Assemble  the  matrix  in  compact  form. 

3.  Choose  a  solution  procedure  and  renumber  the  equations  to  reduce  the  storage 
demands  of  the  factors  of  K. 

4.  expand  the  compacted  K  into  a  data  structure  suitable  for  the  solution  method. 

5.  solve  the  associated  system  of  equations. 

When  there  is  insufficient  primary  storage,  the  assembly  of  the  matrix  in  compacted  form 
opens  a  number  of  possible  avenues  that  one  can  take  to  reduce  the  I/O  cost.  The  expansion  of 
the  compacted  form  need  not  be  done  immediately  after  it’s  assembly.  The  matrix  can  be  kept  in 
compact  form  and  put  on  secondary  store  and  expanded  into  a  full  profile  form  only  when  a  fac¬ 
torisation  must  be  performed.  This  way  the  number  of  data  entries  that  is  read  (in  the  Input 
phase  of  I/O)  can  be  reduced  considerably,  which  in  turn  results  in  a  reduction  in  the  solution 


time. 


Figure  1. 


Finite  element  mesh  of  example  1. 


Figure  2.  Nodal  graph  for  the  mesh  of  example  1. 
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Figure  3.  The  graph  of  unknowns  for  the  mesh  of  example  1. 
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Example  2:  Cantilever  Struc¬ 
ture,  left  end  fixed,  plane  strew 
elements  with  2  degrees  of  free¬ 
dom  per  node. 


Example  3:  Small  Cylinder, 
both  ends  fixed  in  tangential 
direction,  plane  strain  elements 
with  2  degrees  of  freedom  per 
node. 


Example  f  Large  Cylinder, 
both  ends  fixed  in  tangential 
direction,  plane  strain  elements 
with  2  degrees  of  freedom  per 

nc.de. 


Example  5:  4X4X4  Solid  cube, 
fixed  base,  solid  elements  with  3 
degrees  of  freedom  per  node. 


Example  6:  8X8X8  Solid  cube, 
fixed  base,  solid  elements  with  3 
degrees  of  freedom  per  node. 


•SSH 


Table  6.  Description  of  test  ex  nm pies. 
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Appendix  At  Program  Listing 

SUBROUTINE  ELCNT( NUMNP , NUMEL , NEN , NEN 1 , IX. 1C) 

DIMENSION  I X( NEN1 , 1 ) , I C( 1 ) 

INPUT 

NUtffJP  TOTAL  NO  OF  NODES  IN  THE  MESH 

NUMEL  TOTAL  NO.  OF  ELEMENTS  IN  THE  MESH 

NEN  MAX  NO  OF  NODES  PER  ELEMENT 

NENI  DIMENSION  OF  IX  ARRAY 

IX  ELEMENT  CONNECTIVITY  ARRAY 

OUTPUT 

1C  ARRAY  OF  LENGTH  NUMNP.  IT  FIRST  HOLDS  THE  ELEMENT  DEGREE 

OF  EACH  NODE,  THEN  BECOMES  A  POINTER  FOR  AN  ARRAY  THAT 
CONTAINS  THE  SET  OF  ELEMENTS  CONNECTED  TO  EACH  NODE. 

COUNT  THE  NUMBER  OF  ELEMENTS  EACH  NODE  BELONGS  TO 

CALL  IZERO( IC.NUMNP) 

DO  1  10  N  =  1  .NUMEL 
DO  100  J  =  1 , NEN 
I  =  IX( J ,N) 

!F( I  GT  0)  IC( I )  =  IC( I )  +  1 
CONTINUE 
CONTINUE 

SET  UP  POINTERS 

DO  1  20  I  =  2  .NUMNP 

IC( I )  =  IC( I )  +  IC( 1-1 ) 

CONTINUE 

RETURN 

END 


SUBROUTINE  CASSEM( D , A , B , S , P , JCOLE , IRCW.LD, ID .NST.NEL , AFL , BFL ) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

LOGICAL  AFL. BFL 

D I MENS  I ON  D(  1  ) , A( 1 )  ,  B(  1  ) , S(NST, 1 ) ,P( 1 ) . JCOLE( 1 ) . IRCW(  1 ) ,LD( 1 ) 

»  .ID{1) 

COMPACT  ASSEMBLY  OF  PROFILE  MATRIX 

DO  200  J  =  1 .NEL 
N  =  LD( J ) 

IF  (  AFL  AND  N  GT  1  )  THEN 
DO  150  1  =  I , NEL 
K  =  LD( I ) 

IF  (  K  GT.  0  AND.  K  .  LT.  N  )  THEN 

INZ  -  !NZA(  JCOLE (N-  I  )+ 1  ,  JCOLE(N)  ,  IRCW.K) 

A ( INZ)  =  A( INZ)  +  S(  I  ,  J  ) 

END  IF 
CONTINUE 
END  IF 

ASSEMBLE  THE  DIAGONAL 
IF  (  N  GE  I  )  THEN 

IF  (  AFL  )  D(N)  -  D(N)  +  S(J,J) 

ASSEMBLE  THE  LOAD  IF  NECESSARY 
IF  (  BFL  )  B(N)  =»  B(N)  ♦  P(J) 

END  I F  a 

CONTINUE 
RETURN 
END 


o  ci  o  ocio^n  *»  ►*  (innnoooooooo 
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SUBROUTINE  CCMPRO(NUk»JP  ,  NUMEL  ,  NEN ,  NEN1  ,  NDF  ,  IX,  ID,  1C,  IROW,  IELC. 
l  JCOLE ,KP) 

DIMENSION  I  X ( NEN 1 , 1 ) . 1D(NDF, 1 ) , IC( 1 ) , IRCW(  1 ) , !ELC( 1 ) . JCOLE( 1 ) 

FOR  (NUWP, NUMEL ,NEN,NENl , IX, 1C)  SEE  SUBROUTINE  ELCNT 
. . . .  INPUT 

NDF  NUMBER  OF  UNKNOWNS  AT  EACH  NODE 

ID  ACTIVE  UNKNOWNS  AT  EACH  NODE 

. , . .  OUTPUT 

I ELC  HOLDS  THE  SET  OF  ELEMENTS  CONNECTED  TO  EACH  NODE 
1ROW  ROW  NUMBER  OF  EACH  NONZERO  IN  THE  STIFFNESS  MATRIX 
JCOLE  END  OF  ENTRIES  IN  I  ROW  FORM  A  GIVEN  COLUM4 

FIND  ELEMENTS  CONNECTED  TO  NODES 

CALL  I  ZERO  ( IELC, IC(NULNP) ) 

DO  no  N  =  I  .NUMEL 
DO  220  J  =  1 , NEN 
I  =  IX(J.N) 

IF  (  I  GT  0  )  THEN 
KP  =  IC(  I  ) 

00  IF  (  IELC(KP)  EQ  0  )  GO  TO  210 

KP  =  KP  -  I 
GO  TO  200 

10  1ELC(KP)  =  N 

END  IF 

20  CONTINUE 

30  CONTINUE 

SET  UP  COMPRESSED  PROFILE  POINTERS 

KP  =  0 
NEP  =  1 

DO  350  I  =  1  ,  NUMNP 
NE  =»  I C  (  1  ) 

DO  340  1  1  =  1  , NDF 
NEQ  =  ID(  I  I  .  I  ) 

IF  (  NEQ  GT.  0  )  THEN 
JCOLE (NEQ)  =  KP 
KFO  =  KP  +  1 
IF  (  NEP  LE  NE  )  THEN 
DO  330  N  =  NEP ,NE 
NN  =  IELC(N) 

DO  320  J  =  I , NEN 
K  =  I X ( J ,NN) 

DO  3  10  J J  =  1  , NDF 
NEQJ  -  I D ( J  J , K ) 

IF  (NEQJ  GE  NEQ  OR  NEQJ  LT  0)  GO  TO  310 

CHECK  TO  SEE  IF  NODE  ALREADY  IN  LIST 

IF  (  KPO  LE  KP  )  THEN 
DO  100  KK  =  KPO, KP 

I F (  IROW(KK)  EQ  NEQJ  )  GO  TO  310 
300  CONTINUE 

END  IF 
KP  =»  KP  -*•  I 
IROW(KP)  -  NEQJ 

310  CONTINUE 

320  CONTINUE 

330  CONTINUE 

JCOLE  (NEQ)  =*  KP 
END  IF 
END  IF 

300  CONTINUE 

NEP  *»  NE  +  I 
350  CONTINUE 
RETURN 
END 
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INTEGER  FUNCTION  1NZA(N1 ,N2 , IROR.K) 
DIMENSION  IRCM>*(I) 

FIND  THE  TERM  FOR  THE  ASSEMBLY 

DO  100  N  =  Nl ,N2 

IF  (  !ROW(N)  EQ  K  )  THEN 
INZA  =  N 
RETURN 
END  IF 
100  CONTINUE 
C  ERROR  IF  LOOP  EXITS 

STOP 
END 


SUBROUTINE  IZERO(IA.NN) 
DIMENSION  IA(NN) 

DO  100  N  =  1 ,NN 
I  A(N)  =  0 
100  CONTINUE 
RETURN 
END 
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